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Abstract 
The  structure  cf  the  interface  of  an  argcn-like  fluid  in 
ecuilibriiT  with  its  own  vapor  st  low  terrperature  is  studied 
using  rolecul=>r  dynarrics.   The  lonritudinal  pair  correlations  in 
the  interface  are  found  to  be  consistent  with  a  sirply  defined 
enserble  of  local  therrodynaric  states.   However,  the  transverse 
correlations  exhibit  very  long  range  behavior  not  predicted  by 
straightforward  local  therr^cdynarics .   These  results  strongly 
suggest  that  the  interface  is  rade  up  cf  an  enseirble  of  configura- 
tions in  each  of  which  the  transition  froir  liouid  to  vapor  is 
locally  sharp,  but  that  the  transition  surface  fluctuates  strong- 
ly in  space  and  tire. 


1 .   Introduction 
Since  the  inception  of  statistical  mechanics  the  structure 
of  the  transition  zone  between  liquid  and  vapor  in  equilibriun 
has  been  of  great  theoretical  interest.    From  a  thermodynamic 
viewpoint,  surface-associated  Quantities  such  as  surface 

tension  fit  without  difficulty  into  a  framework  encompassing 

2 

the  variation  of  both  bulk  and  surface  parameters.    But  the 

jump  from  a  microscopic  description  of  structure  in  bulk  fluid 
to  that  in  the  interface  repion  has  long  eluded  reliable 
analysis.   This  is  largely  because  the  experimental  information 
so  necessary  to  direct  theoretical  efforts  remains  very  meager. 
There  are  relatively  few  eouivalent  particles  entering  into  any 
given  aspect  of  the  interface,  and  their  contributions  to 
physical  measurements  are  marked  by  the  background  from  bulk 
particles . 

The  advent  of  modern  digital  computers  has  changed  this 
situation.   Now  it  is  possible  to  perform  a  numerical  simulation 
of  an  equilibrium  system  with  a  large  proportion  of  interface 
particles,-  and  sufficiently  many  to  ensure  good  statistics 
not  only  for  density  computations,  but  even  for  joint  distribu- 
tions of  particle  pairs.   With  this  information,  one  can  begin 
to  evaluate  various  theoretical  approximations  and  suggestive 
concepts,  and  more  importantly  provide  a  suitable  qualitative 
framework  for  the  microscopic  phenomenology  of  two-phase  fluid 
interfaces. 
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In  this  paper,  v;e  will  study  an  argon-like  fluid  at 
low  temperature.  It  is  now  fairly  well  established   that  for  a 
planar  two-phase  interface,  the  density  profile  p(z)  is  a  smooth 
monotone  function  approaching  the  corresponding  liauid  or  vapor 
density  in  the  one-phase  regions  (see  Fig.  1  for  profile  from 
data  of  this  paper).  The  common  ingredient  of  theoretical 
derivations  of  p(z)  has  always  been  some  ansatz  for  relating  the 
pair  distribution  function  p  ^  (r-,r„)   to  p(z)  and  a  suitable 
bulk  pair  distribution.   In  the  present  analysis,  we  first 
investigate  the  longitudinal  pair  correlations  and  show  that 
they  are  consistent  with  a  simply  defined  ensemble  of  local 
thermodynamic  states.   Self-consistent  use  of  this  relation 
however  produces  a  density  profile  substantially  steeper  than 
the  numerically  observed  p(z). 

V/e  next  proceed  to  transverse  correlations,  via  a  trans- 
verse Fourier  transform,  and  find  that  they  exhibit  very  long- 
range  behavior  not  predicted  by  naive  local  thermodynamics. 
This  behavior  is  shown  to  be  a  ouite  general  consecuences  of 
the  existence  of  a  self-supported  phase  transition  region,  and 
is  placed  dynamically  in  the  context  of  ensembles  of  wave-like 
gas-liquid  interfaces.   The  observed  broadening  of  the  density 
profile  has  a  similar  origin.   To  exhibit  the  genesis  of  the 
transverse  correlations,  an  intermediate  energy-density  contour 
of   each  configuration   of  the  computer-generated  system 
is  carried  out.   The  results  strongly  reinforce  the  picture  of 
an  ensemble  of  configurations,  each  consisting  of  sharply 
divided  liquid  and  gas  regions,  but  with  a  highly  fluctuating 
interfacial  surface. 


2.   Longitudinal  Correlations 
Our  model  system  consists  of  1728  particles  interacting 
via  a  classical-rriechanical  Lennard-Jones  (6,1?)  potential 

V(r)  -  V(rp)  ,        0  <  r  <  rp 
=  0  ,  r^  <  r  (1) 

where   V(r)  =  '^^{{^/r)'^^   -    (c/r)^) 


Me   choose  ^ ,  ^,    and  rr  =  (r.^/^Sf.)      '^   as  length,  energy,  and 
tire  units  (for  real  argon,    „  =  3-^05  A  ,    e  =  119.'4  K, 
Tp  =  3.112  X  1C~  ^   sec).   The  particles  are  placed  in  a  fully 
periodic  box  of  size  L  x  L  x  3L  (with  L  =  13.15a),  with  a 
Maxwellian  distribution  of  velocities  corresponding  to  a  mean 
temperature  of  0.70^  (8^°  K).    The  details  of  creating  an 
inhorrogeneous  systerr  without  any  external  field  in  this  periodic 

c 

box,  using  rrolecular  dynairics,  are  presented  elsewhere.'' 

After  the  filin  eauilibrates  (showing  fairly  constant  temp- 
erature and  pressure  profiles),  the  configurations  are  stored  on 
rr-apnetic  tape   at  intervals  of    0.128xp,  and  equilibrium  quanti- 
ties computed  as  averages  over  800  configurations,  as  well  as 
over  any  theoretical  directions  of  invariance.   The  full  density 
profile  that  develops  is  shown  in  Figure  2,  consisting  of  a 


vapor  of  bulk  density  n-  =  .0033'   Here,  interfaces  are  perpendi- 
cular to  the  long  (z)  direction. 

The  shape  of  the  density  profile  has  been  the  subject  of 
numerous  theoretical  studies,  mainly  before  reliable  numerical 
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knowledge  existed.   Each  of  these  was  corrpelled  to  nake  a  critical 
assumption  as  to  the  forip  of  the  pair  distribution  function  in 
the  inhornogeneous  system.  We  are  now  in  a  position  to  evaluate 
the  assumptions.   Since  the  mean  density  variation  is  only  in  the 
longitudinal  (z)  direction,  we  first  concentrate  upon  the  longitudi- 
nal correlations.     The  two-particle  distribution    p    (Ci'C?) 
denotes  the  number  of  distinct  pairs  of  particles  per  unit  volume 
in  r^-space  and  unit  volume  in  rp-space.   But  our  system  is 
constructed  to  be  translation-invariant  in  the  x  and  y  direc- 
tions, so  we  can  write 

p^^^r^,!;^)  =  p^^^z^,Z2;x^-X2,y^-yp),  (?) 

2  2  1/2 

indeed   depending   only   on    [(x.-Xp)   +  (y.-y-)  ]      for 

infinite  interface  cross-section.   By  the  loneritudinal  pair  distri- 
bution, we  now  mean 

p^^2\z^,Z2)  =  p^^^z^,z„;0,C).  (?) 

The  assumptions  made  as  to  the  structure  of  an  inhornogeneous 
fluid  of  course  tend  to  follow  local  thermodynamics,  mimicking  -- 
in  the  small  --  successful  conceptual  aspects  of  bulk  systems. 
Microscopically,  the  leading  level  of  structure  for  a  homogeneous 
component  of  density  n  is  given  by  the  radial  distribution  function 

g(r^-r2;n)  e  p  ^  ^^  r  ^-r^  ;n  ) /n"  .  (4) 

How  then  is  this  information  to  be  used  to  estimate  the  inhomogeneous 
P^^^(r^,r2)? 

Consider  for  reference  a  typical  two-phase  PV  isotherm 
-H- 


(Figure  3).   Suppose  that  an  isolated  planar  two-phase  inter- 
face is  defined  --  with  sorre  suitable  recipe  --  by  gas  for 

z  <  0,  liauid  for  z  >  0.   Then  in  the  zeroth  order  ansatz  of 

6 
Nazarian,   one  assumes 

where   n 


<  0 
>  0 


(5) 


effectively  restricting  consideration  to  the  terrrinal  pure  phase 

points  A  in  Figure  3,  the  choice  dictated  by  the  center  of  mass 

7 
of  r.  and  r„.   Toxvaerd   avoids  the  discontinuity  in  n  by 

choosing   n   as  the  center  of  mass  density   and  computing 

g(r.-rp;n)  by  extrapolation  along  the  metastable-unstable 

Van  der  Waals  loop  BE: 


(2) 


(r  ,r^)  =  p (Ci )p (Cp^g^Ei-CoJp fi  (z  +z^)]) 


(6) 


To  the  extent  that  (6)  is  a  valid  starting  point,  first  order 
density  gradient  corrections  can  be  shown  to  vanish.    A  third 
alternative,  which  we  shall  examine  in  some  detail,  is  to  move 
along  the  two-phase  coexistence  line  C.   This  requires  some 
explanation . 

For  the  bulk  fluid  in  two-phase  coexistence,  statistical 

g 
mechanics  as  ordinarily  practiced   attaches  a  probability 

Pq  =  1  -  p^  at  each  space  point  for  the  system  to  be  liquid, 

and  computes  all  distributions  with  this  weight: 

n  =  pp  no  +  p^n^  , 
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The  irplication  is  that  gas  and  liquid  are  predoipinantly  present 
only  in  very  large  clusters,  so  that  in  the  course  of  time 
a  given  point  --  or  set  of  points  --  can  be  regarded  as  either 
internal  to  a  liquid  cluster  or  internal  to  a  gas  cluster. 
If  attention  is  focused  on  spatial   and  teirporal  regions  of 
transient  existence  of  small  clusters,  then  the  interfacial 
free  energy  affects  the  local  pressure  which  can  be  supported, 
and  the  non-negligible  interfacial  volume  must  be  averaged  into 

the  environment  of  a  given  point.   This  is  presumably  what  (6) 

1 0 
takes  advantage  of.   But  there  is  evidence    that  the  struc- 
tural detail  in  the  vicinity  of  a  rapid  density  change  is  very 
short  range  and  in  particular  on  a  scale  small  compared  to  that 
of  Figure  1.   Hence  a  suitable  dynamic  picture  of  a  two-phase 
system  should  indeed  include  the  possibility  of  an  interface 
sweeping  past  a  given  spatial  point,  but  this  point  can  still 
be  regarded  at  any  instant  as  either  in  gas  or  in  liquid. 

The  interpretation  of  the  foregoing  is  clear  enough. 
We  can  view  the  two-phase  system  as  a  temporal  sequence,  or 
phase  space  ensemble,  of  configurations  each  sharply  divided 
into  gas  and  liquid  regions.   The  geometric  form  of  the  dividing 
surface  is  however  open.   As  a  primitive  approximation,  and 
in  accordance  with  the  resulting  symmetry,  we  shall  take  this 
as  planar.   For  a  given  division  plane,  positioned  at  Z,  we 
have  pure  gas  for   z  <  Z,    and  uncorrelated  pure  liquid  for 
2  >  Z.  It  follows  that 
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P^Cr)  =  Dq  e(z  -  Z)  +  n^e(Z  -  z) 

P^'^'^r^.r^)  =  n/g^Cr^-r^)  t  (z  ^-Z  )  e  (z.-Z  )       ^^^ 
+  npn^(1-e(z^-Z)e(z^-Z)  -  e(Z-z^)  eCZ-z^)) 


n^^^Si  ^{r  ^-r^)    e(Z-z^)  cCZ-z^)  , 


where  £(y)  is  the  Heaviside  step  function,  i,e.    (1  +  sgr(x))/2. 

Now  if  Z  is  distributed  with  probability  density  f(Z)  and  corres- 

f  Z 
ponding-  cumulatiT/-e  F(Z)  =  I  ffZ')  dZ',   we  have  on  integration 


p^^^r^.r^)    =    Hp^gpCr^-r^)    F(Hin  (z  ^  ,  z„  )  )  (9) 

+   n^n^CFCMaxCz^.z^))    -   F(Min    (z^,z^)) 
+   n^^g^(r^-rp)(1    -   F(Hax    (z^.z^)). 

( 2 ) 
The  desired  relation  between  p    (r.r-)  and  p(r) 

results  fror  (o)  on  eliminating  the  function  F(z): 

p^'"  {r  yV^)    -   np  /(n^-n^ )(  P  (Min(z^  ,Zp)  )-n^  jgp(r  ^-r^) 


+  npn^/(np-n^){p(Max(z^,Z2))-p(Min(z^,Z2))J 
+  n^^'/Cnp-n^ )  (np-p(Max(z^  .z^)  )  )g^(r^-r2) , 

(10) 
quite  different  in  forp  frorr  (5)  and  (6).   However,  if  n^  can  be 
neglected,  as  is  the  case  with  our  data,  (10)  reduces  to 

o^^hvyV^)    =  npP(Min(z^,Z2))  g^{v^-r^))  (10') 
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very  ff'uch  in  the  same  mold  as  (5)  and  (6).   The  radial  distribu- 
tion function  at  liquid  density  to  serve  as  input  to  (10')  is 
shown  in  Figure  4. 

We  can  now  return  to  the  question  of  longitudinal  correla- 
tions.  According  to  (10')  we  should  have 

PL^2\z^,Z2)  =  npp^.^  gp(Az)  ,  (n) 

(?) 

and  so  the  ratio  p^        ''"o^tnin  ^^^   heep    plotted,  in 

Figure  5,  for  several  values  of  az  h  z^  -  Zp  ,  as  a  function 
of  z  =  (z.  +  z-)/2.   Although  the  corputer  values  of  p.    vary  by 
a  factor  of  50  in  the  range  examined,  this  ratio  remains  very 
nearly  constant  --  and  with  correct  constant  --  at  each  fixed  value 
of  AZ. 

We  conclude  that  the  approximation  represented  by  an  ensemble 
average  of  sharply  separated  gas-licuid  regions,  even  with  the 
further  assumption  of  a  strictly  planar  interface,  provides  an 
adequate  picture  of  longitudinal  correlations  in  our  model  system 
with  its  model  parameters.   A  similar  test  of  (6)  requires  much 
more  bulk  system  input  data  and  is  correspondingly  less  definitive, 
but  does  not  lead  to  substantially  poorer  results  in  the  region 
examined.   The  advantage  of  (10)  and  (10')  is  their  greater 
simplicity  and  deeper  dynamical  implications.   In  addition,  they 
can  be  easily  improved.   We  will  pursue  these  advantages. 
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3.   The  Density  Profile 

(?) 
The  conventional  way  of  using  inforrration  on  p    to  obtain 

an  interphase  density  profile  has  been  to  insert  this  into  an 

exact  relation  connecting  o^    and  o.   For  this  purpose,  the 

first  of  the  EEGKY  hierarchy  (local  balance  of  mechanical  and 

thermodynaDic  forces) 

7p(r)  +    e  I  p^^^(r,r')  V4(r  -  rM  d-r'  =  0^       (12) 


with  ^   the  pair  interaction,  has  generally  —  but  not  univer- 
sally  —  been  employed.   The  Qualitative  nature  of  the 
profiles  derived  fror  (H)    or  (6)  is  easily  established.   We 
consider  any  region  in  which  p(r)  is  close  to  the  bulk  fluid  (gas 
or  liquid)  density  n,  and  hence  set  p  '^  (r,r')  =  p(r )  p(r '  )g(r-r '  ) , 
where  g  is  referred  to  density  n.   For  a  plane  syroiretric  profile 
o(z),  (12)  becorres 

p'(z)  =  -Bp(z)  I  g(r' )(z'/r' )*'(r' )p(z-z' )  d'r'.    (13) 

AssuR'ing  g*'  short  range  on  the  scale  of  p,  we  can  Taylor 
expand  p(z  -  z')  about  z,  so  that  on  using  the  non-vanishing 
angular  averages  <z'^^>  =  r'^^/(2s  +  1),  we  have 

P'(z)  =  -8p<z)  j  g(r')  ♦'(r')[-  ^   r'p'(z) 

-  -^   r'^P"  '(z)  +  ...]  d-r'.  (1i4) 

Hence,  inserting  the  virial  pressure, 

p'(z)  =  ap'(z)  +  bP(z)p'"(z)  +  ...  (14') 

where    a  =  2(n-6p)/n^,   b  =  (6/30)  I  g(r • )r • * • (r  •  )r '^  d^r'. 
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Tc  solve  (14'),  we  need  only  integrate,  ottaining 


(z)  =  K  +  0.5  ap(z)'  +  0.5  b(2p(z)p"(z)  -  p'Cz)^)    (15) 


p 
for  suitable  K.   Now  multiply  by   p'(z)/p(z)    and  integrate, 

yielding 

K^  +  In  p(z)  =  -  K/p(z)  +  0.5  ap(z)  +  0.5  b  p'(z)^/p(z)    (15') 
for  suitable  K'.   If  p(z)  =  n  is  to  satisfy  (15)  and  (15'), 
K  and  K'  are  determined,  allowing  (15')  to  be  written  as 


n^b(p'(z)/n)^  =  2p(z)[ln  p(z)/n  -  1  +  n/p(z)] 


Finally,  expanding  in  small  a(z)  h  p(z)/n 


(16) 


n^b(A'(z))^  =  bPa(z)^  .  (16') 


Working  in  from  gas  and  liouid  ends,  this  results  in  exponen- 
tial deviations  into  the  interface  region,  with  a  much  more 
rapid  relative  rise  from  the  low  n  gas  side  --  see  Figure  6. 
It  is  then  the  nature  of  the  connection  region  that  distinguishes 
e.g.  {5)  from  (6). 

The  need  for  simultaneous  knowledge  of  $  and  g  is  more  than 

a  small  inconvenience  if  one  wants  ultimately  to  make  use  of 

(2) 
experimental  data.   However,  one  can  always  write  down  a  p    -  p 

relation  in  which  the  internal  potential  *  does  not  appear  at 

all.   This  is  most  easily  derived  by  putting  the  sytem  in  an 

1  2 
external  field  u(r).   It  is  known    that  at  fixed  chemical 

potential,  a  change  in   u   produces  a  change  in  p  according  to 

the  linear  response  relation 
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6p(r)  =  -B  I  [p^^^(r,r')  +  p(r)6(r-r')  -  p (r) p (r* )  ]«u(r • )d-r • 
i  -B  j  S^^^(r,r')  «u(r')  d'r'  (17) 

or  its  inverse 

-  8«u(r)  =  j  C^^^(r,r')  6p(r;')  d-r'.  (17') 

Here 

C^^^(r,r')  r   6(r-r')/p(r)  -  c(r,r')  ,         (18) 

where  c(r,r')  is  the  usual  direct  correlation  function.   For  a 
translation-invariant  internal  potential,  a  unit  translation 
cf  the  whole  systeir-  produces 

su(r)  =  7u(r),   Sp(r)  =  7p(r)  ,   whence 

-  6vu(r)  =  i  C^2^(r,r')  vp(r')  d^r'.  (19) 

1^  15 
Now  if  u(r)  =  0,  one  has  the  desired   ' 

7p(r)  =  p(r)  j  c(r,r')  Vp(r')  d-r'.  (20) 

Approximations  inserted  into  (20)  rray  of  course  give 
dramatically  different  results  than  when  inserted  into  (12). 
But  the  rrost  convenient  approxiirations  are  in  fact  those  for 
c  ^^  and  not  for  p    .   The  analogs  of  (iJ)  or  (6)  have  not 
yet  been  used  for  this  purpose,  and  indeed  the  function  c  would 
be  expected  to  become   singular  at  infinite  compressibility  in 
the  metastable  region.   However  the  wings  of  the  interphase 
region  can  presumably  be  analyzed  as  was  (12).   We  expand  the 
plane  syrrmetric  case,  obtaining  for  c(r,r')  =  c(r-r') 
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(z)  =  p(z)  j  c(r')  p'(z-z')  d^r' 
=  p(z)  (Ap'(z)  +  Ep"'(z)) 


c(r')  d^r-  =  -^(1  -   -f^^) 

n        an 


(21) 


B  =  I  I  r'^  c(r')  d^r', 
and  conclude  as  before  that 

2nB(p'(z))2  =  (1  +  ^^)  a(z)2  ,  (21') 

with  sirrilar  Qualitative  interpretation. 

(2) 
There  is  of  course  no  point  to  usinp  a  p     -  p  relation 

rrore  accurate  than  the  approximation  to  p     ,    IL   the  sare 

mechanisin   is  responsible  for  the  inaccuracy  of  each.   The 

approximation  (10)  was  based  upon  neglect  of  the  gas-liquid 

interaction  and  the  associated  longitudinal  fine  structure, 

so  that  for  each  configuration  of  the  relevant  ensemble,  each 

volum.e  element  either  belongs  to  liquid  at  density  Hq  or  gas 

at  density  n^.   Suppose  that  we  have  some  microscopic  criterion 


G(r)  >  0    implies    liquid  at  r 

(22) 
Q(r)  <  0    implies    gas    at  r  . 

Then  the  mean  density  at  r  under  the  above  assumption 

becomes  simply 

p(r)  =  nn<e(Q(r))>  +  ni<e(-C(r))> 

°  ^  (23) 

=  0.5  (n^  +  n^)  +  0.5  (n^  -  n ^ )<sgn(0( r ) ) >  . 
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K'ow 


<sg 


n  C(r)>  =  -4   P  I  <e^^°^r)>  ds/s 

ITl        J 


2  2, 


4   P  I  exp  [isC(r)  -    s'^a  (r)lds/s 
1     J 


(2J4) 


to  second  cuniulant  order,  where 


(24') 


Evaluating    (by   applying    3/30),    vje    have 


(r)    =   n 


Q(r)/a(r) 
I       e    ^    ^^    dt 


/27r 


(25) 


For  plane  symrretry,  p,  Q,  o  ^re    functions  of  z  alone.   Suppose 

that  the  mean  interface  is  at  z  =  0,  i.e.  C(0)  =  0.   Then,  to 

leading  order  in  z,  we  have  the  error  function  representation 

z/x 

p(z)  =  n.  +  ^ ^    I  e  ^  ^^  dt  (26) 

^     /2w      J 

where   x  =  o(0)/C'(0). 

But  what  shall  we  choose  as  Q(r)?   In  order  to  make  use  of  one 

and  two-body  data,   C  must  contain  only  0  and  1-body  terms,  and 

of  course  be  invariant  to  translation  of  z  and  of  the  whole  system. 

Thus 

0(r)  =  Zi  T(r.  -  r)  -  tq  .  (27) 


Since  C^(0)  =  0,  we  have  t 


(z* )T(r' )  d-r' ,  and 


so  the  parameters  of  (24')  become 
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P(z)  =  I  [p(z' 


z) 


p(z'))T(r' )  d^r 


a^(z)  =  j  j  S^2)(^,^^„)^(^,_2)^(r'"-z)  d^r'  d-r". 


In  particular 

C'(0)  =  0 
2 


(28) 


(28' 


(C)  =  I  I  S"(r' ,r")T(r')T(r")  d^r'  d'r" 


In  the  nurrierical  case  under  discussion,  the  pas  density  n-  is 
very  low.   Let  us  choose  n.  =  0,  so  that  (10')  applies.   After 
sorre  manipulation,  the  distance  parameter  x  of  (26)  then  takes 
the  form 


,^  =  [  i  i  p(z^.^)SQ^^^r-r')T(r)T(r')/np  r^^r  d^r' 
+  I  I  ^"o-P^^max^^''^^min^^^C)^(C')]/U  p'(z)T(r)  d -V )  ^ , 


(29) 


and  our  task  is  to  solve  this  for  x 


Issuming  that  the  range  of  S 


(2) 


is  srall  on  the  scale 


of  interest,  we  shall  make  the  replacement   S^(r-r') 


I  s 


(2) 


(r")  d^r"  6(r-r') 


anp/aeP   in  equation  (29) 


Finally,  there  is  the  Question  of  the  choice  of  the  density  test 
function  t,  although  the  results  should  not  be  sensitive  to 
this  choice.   We  shall  take  t  as  Gaussian: 

1        -  I   r'/u' 


(r)  =  T(x)T(y)T(z)  = 


(30) 


so  that  Q(r)  is  indeed  a  measure  of  the  local  density  at  r. 
Now  (29)  can  be  written  as 


x^(  I  p'(z)  t(z)  dz)^  =  <n„  I  p(z)t(z)^  dz  (j  T(z)2dz)' 


(31) 


I  I  p(z)p(z')t(z)t(z')  dz  dz 
;  J 


z+z  '<0 
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a  somewhat  tedious  evaluation,  (31)  achieves  the  forn' 

X-  =    ^-^^  <  /(2F(x)) 

wliere   F(x)  =  x^  [4,^/ (  1+x^)-x^(ir/i4  +  tan"^  x//2+x^)] 
and     x  =  v/x  .  (32) 

To  rrake  sure  that  x  is  insensitive  to  the  value  of  v,  we  demand 
that  F'(x)  =  0,  resulting  in  the  numerical  value  F(x)  =  3^  at  x 
=  v/x  =  1.6.   Hence 

1/3 

(33) 


X  -^  0.  4^4  < 


(0)  '^    n^J/T^    X 


0.92  n 


-1/; 


This  approximation  should  be  appropriate  to  the  center  of  the 
interface  region  rather  than  the  wings.  Using  the  measured  value 
<    -    0.04,  we  find  from  (33),   p'(0)  '^    2.1.   On  the  other  hand, 
the  measured  density  profile  is  much  broader,  with  p'(0)  '■^    0.35. 
We  shall  soon  see  why  this  is  the  case. 

4.  Transverse  Correlations 
Neither  version  (6)  nor  (10)  of  local  thermodynamics  appear 
to  show  major  qualitative  deficiencies  when  longitudinal  effects, 
be  they  density  profile  or  density  correlations,  are  in  question. 
The  situation  deteriorates  rapidly  but  informatively  when  trans- 
verse correlations  are  considered.  For  pure  transverse  correla- 
tions, we  identify  the  z-values  of  the  points  in  question,  and 
thus  define  the  transverse  pair  distribution  by 
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(2) 
where 


(2; 


(x^-x^,  y^-y2) 


(34) 


According  to  approximation  (6) 


>^^^z,r^2)  "  P^^)^  ^^Ci2'  ''^^^^ 


while  accordinf  to  (10) 


(2) 


P (z)-n 


T   ^^ 'Ci2^ 


-  "ngn^r'io) 


7—^ p,F,(r,^)  . 


(35) 


(36) 


The  corresponding  aoproxiir'aticns  for  the  transverse  correlation 


(2) 
function    F^    ^(z,r^2) 


(2) 


(z,r.p)    -    p(z)      then    becorrie 


T 
,(2) 


F^^^r,,;p(z)) 


12 
p(z)-n 


F^2)(^       ) 


'o~°^^^     (21 

-L F       M  r       ) 


(37) 


(38) 


It  is  convenient  to  go  over  to  transverse  Fourier  space: 

ik.r , ^   ^ 


p(z) 


I  F 


(2) 


(z.r^^)  e 


12 


(39) 


While  (37)  and  (38)  evoke  similar  behavior  in  F^(z,k)  for  k  /  0, 
i.e.  a  bulk  F(k)  at  some  mean  density  compared  to  a  mean  of  Fg^if^ 
and  F.(k),  the  behavior  at  k  =  0  is  qualitatively  different: 
regular  for  (37)  but  singular  for  (38).   How  do  these  compare  with 
our  numerical  simulation?   In  the  simulation,  we  compute 
N(Az)  -i^'(Zl~Z2^ 


z 


>/<N(Az)> 


(40) 


where   k  =  (2Trn/L,  2Trm/L,0),   m,n  =  0,1,2,...;  mn  ^    0, 
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i  and  j  running  over  the  NC^z)  atonis  in  a  slice  of  volutre 

L  X  L  X  AZ   centered  at  z.   The  angular  brackets  denote  the 

2 
average  over  all  configurations,  and  corrrDon  values  of  k  are 

averaged  as  well.   For  sirall  iZ,  it  is  easy  to  see  that 


for  two  values  of  z,  respectively   in  the  rr,iddle  of  the  bulk 
liquid  and  in  the  middle  of  the  interface.   In  Figure  8,  the 
circles  show  S„(k)  calculated  with  1024  particles  in  a  box 
of  size  6.78a  ^  6.78a  x  56.2^40  .   The  close  agreement  down  to 
the  smallest  k  of  the  smaller  system  lends  credence  to  the 
reliability  of  the  numerical  results. 

The  enhancement  of  low  k  values  in  the  interface  region 
is  brought  out  strikingly  in  Figure  8.   There  is  apparent 
divergence  as  k  -►  0,  indicating  transverse  correlations  over  the 
size  of  the  box.   One  can  say  that  the  totally  correlated  k  =  0 
component  of  (38)  is  in  reality  smeared  into  a  merely  very  long 
range  correlation.   To  bring  out  the  contrast,  we  plot  in  Figure 
9  the  difference  between  S_  obtained  from  simulation  and 
that  predicted  by  (38),  for  the  minimum  value  of  k.   This 
difference  is  indeed  confined  to  the  interface  region,  and 
closely  resembles  the  curve  of  p'(z)'^,  also  drawn.   We 
will  comment  on  this  relation  in  a  moment.   Lack  of  reliable 
knowledge  of  the  predictions  of  (37)  preclude  a  similar  compar- 
ison.  But  we  can  say  in  general  that  its  extrapolation  from 
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physical  (one-phase)  densities  will  certainly  not  introduce  the 
observed  singular  k  behavior.   On  the  other  hand,  an  independent 
assessment  of  the  rretastable  state  correlations  should  at  least 
exhibit  a  divergence  as  k  +  0  at  the  points  of  infinite  corr^press- 
ibility. 

The  genesis  of  singular  low-k  transverse  behavior  has 


been  suggested  in  a  recent  paper  of  Wertheim 


15 


A  brief  step 


beyond  his  analysis  is  all  that  we  require.   Suppose  that  we 

(2)       (2) 
take  the  transverse  Fourier  transforrr'S  of  C    and  S    , 

denoting  these  by  C  and  S.   Then  (see  (17,17'), 


I  C(z.,z_;k)  S(z2,z,;k)  dZp  =  6(z.-Zp), 


(ill) 


so  that  C(k)  and  S(k)  are  again  inverse  matrices.   On  the 
other  hand,  according  to  (19,20)  in  the  plane-syrrrDetric  case. 


-(2) 


(z.  ,  z,  ;r 1 ,)  p  •  (z,)  d-r,  =  0  , 


12 


(M2) 


I  C(z^  ,z^;0)    p  '(z^)  dz. 


0  . 


(ii2') 


Thus  C(0)  has  p'(z)  as  a  zero-eieenvalue  eigenfunction ,  and 

2 
does  not  have  an  inverse.   But  if  C(k)  is  a  function  of  k 

2         2 

alone,  analytic  in  k   around  k   =  0 ,  we  can  certainly 

write 


C(z,,z„;k)  =  C(z,,z„)  +  k  C'(z.,z) 


I  C(2)( 


1'"2;':i2)'^12  ^^^^2^ 


(^3) 


In  a  representation  including  the  C(0)  eigenfunction 
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,,,.(z)  =  d'(z)/[  I  (d'(z))^  dz]^/2  (HH) 

^  1 

(assumed  now  degenerate),  the  only  spall  diagonal  eletrent  of  C(k) 
is 

C„n(k)  =  ^^    I  p '(Zi)C'(z.,Z5)  p'(z„)  dz.  dz„ 
/  j  p  '(z)^  dz  . 

Thus  (readily  verified  by  standard  perturbation  theory)  the  only 

singular  element  of  S(k)  in  Spp(k)  =  1/Cpp(k)   and  we  have 

.  D  '  (z  -  )  p  '  (Zn) 

S(z^,Zp;k)  =  —^ ■ +  ... 

k^  I  p'(z)  C'(Zt,z„)  p'(z^)  dz^  dzp 

together  with  non-singular  terrrs  in  k. 

The  p'(z)^/o(z)  dependence  of  F„( 

indeed  that  observed  in  Figure  6.   Only  the  magnitude  of  the 

singularity  remains  to  be  discussed.   Eut  in  fact  the  surface 


^  j  p'(z^)r2„  C^^^z^,z^;r^^)p'(z2)dz^dZ2d2r^2-    ^^^^ 


V/e  conclude,  on  identifying  the  deviation  from  a  plane  inter- 
face ensemble  with  the  singular  part  (il5),  that  to  leading  order 

ink,  P 

T   1   p  '(z) 

F^(z,k)  =  7-  ^  .  (^7) 

^    -  '-'    V~    p(z) 

In  the  present  case,  s  =  1/0.7,  y  =  0.7,  p(z)  =  0.35.   If  the  choice 

p'(z)  =  0.50  is  made,  well  within  the  margin  of  error,  the 

simulated  F^ ,  from  (40')  falls  quite  well  on  the  curve  (4?) 

--  see  Figure  10. 
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EQuation  (iJ5)  was  derived  on  auite  general  grounds  of 
consistency,  but  its  forrr  is  quite  suggestive  of  a  ohysical 
mechanism.   The  point  is  this.   AccorOing  to  (38),  there  is  a  6 
function  singularity  in  F,j(z,k)  at  k  =  0.   The  basis  of 
this  is  the  assumption  of  a  plane  interface  between  liauid  and 
vapor.  This  assumption  of  perfect  coherence  is  most  unphysical 
for  large  separation  in  the  interface,  and  we  expect  some 
weaker  small  k  behavior.   Indeed,  in  the  absence  of  infinite 
range  correlations  (as  in  a  perfect  solid),  one  must  have 
the  asymptotic  results 


(2) 


(^1.^2^  *  P  (z-)  )  0  (^2^ 


F^^^z.r^^)  "  ^ 


(^8) 


as  |r^-rp|  -►  <=  or  r^-  -  »  respectively.   These  are  of 

course  consistent  with  (6).   Thus,  although  the  local  z-profile 

may  be  represented  by  (10),  global  structure  is  not.   But 

the  amplitude  of  the  deviation  goes  as  kT/y,  strong  evidence 

that  one  is  looking  at  modes  thermally  excitea  against  surface 

1 7 
tension.   Surface  modes    are  clearly  the  leading  candidates. 

It  is  tempting  then  to  build  up  the  true  surface  configura- 
tion by  "unfreezing"  surface  modes  on  a  bare  or  intrinsic 
surface  profile.   This  however  requires  firm  knowledge  of 
mode-mode  correlations  and  hence  may  not  be  useful  for  detailed 
questions  about  surface  structure.   If  instead  one  looks  at 
infinitesimal  changes  in  surface  mode  amplitudes,  the  objection 
does  not  hold.  Indeed,  pair  correlations  can  be  assessed  in  this 
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f  ? ) 
fashion.   S    (r,r')   in  a  classical  fluid  serves  not  only  as 

density-density  correlation,  but  also  --  see  (17)  —  as  linear 

response  to  ar  inf initesiiral  applied  -field.   Suppose  then  that 

an  external  potential  u(r)  is  applied,  resulting  in  a 

distortion  of  the  z  =  0  surface  of  amplitude  e:(x)  h  c(x,y). 

In  ether  words,  the  new  density  profile  is 

p  (r)  =  o(z  -  5(x)  ) 

Now  --  at  least  at  long  wave  length  --  a  surface  tension  y 

f  2  1/22 

irrplies  a  surface  energy  of  E  =  ^  I  ( ''  +  l'5|  )    <i    ^• 

Hence  the  energy  change  through  second  order  in  u  and  c  is 

aF  =  (y/2)  j  ]V5|^  d^x  +  j  (  p(z)-F,(x)  p' (z)  )u(r)d-r 

=  (./2)(1/L^)  I  k^  ,^,_^    ^    r  ^(^)^(^),?, 

~      ~         '  -   ,  (50) 

-  {^/L    )     I  I  ?L,  D'(z)u(r)€xp[-ik.x]  d^r. 

J       K  ~  ~   ~ 

Using  the  Eoltzrnann  weight  e~^^^,  the  free  energy  change 
associated  with  (50)  becorr'es 

aF(u)  =  I  p(z)u(r)  6'r 

J      .       r  -u  7      9         ?  ^51) 

-  (1/2L^)  z  I  I  p'(z)u(r)e~^~  ~  d^rr/yk^. 

But  Ap(r|u)  =  6AF(u)/6u(r) ,  and  S(r,r')  =  -  6 Ap ( r  |  u ) / « gu ( r ' ) . 
It  follows  that  at  long  wave  length 


and  hence 


;(r,r)  =  (1/6yL^)  [p'(z)p'(z')  e"^^ '^^"^ ' ^ ]/k^ , 


F^(z,k)  =  p'(z)^/[BYk^D(z)] 


(52) 


(52') 


exactly  reproducing  (^7). 
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When  is  it  reasonable  to  exarine  not  just  inf initesiral 
changes  in  S  but  rather  the  full  amplitude  of  surface  waves 
in  the  bare  interface  alluded  to  above?   Perhaps  if  only  s 
one-particle  property  such  as  the  density  profile  itself  is  in 
question.    Suppose  then  that   p  p,  ( z )  =  npe(z)    denotes  the 
idealized  bare  density  profile,  Y„  the  ccrrespondinp-  surface 
tension,  and  ^(x)   the  full  rricroscopic  surface  distortion. 
We  then  have 


p(z)  =  <P^(r)> 


;e(z  -  C(x) 


(5: 


P'(z)  =  ^  f  e^'^^-^'^V^    ds 


(5^) 


According  tc  (50)  at  u  r  C,  <e 


-1S£ 


:exp[-(is/L^)    z    I,     e"^~*~]> 


=    exp    [-   j       (s'^/bycL^)    2     1/1<    ]     •       Hence 


(z)    = 


exp    [-  - 


where 


1 


z^/x^] 


(55) 
(55-) 


^0 


The  k  =  0  terrr  in  (^5')  cannot  of  course  be  present  in  the  forir 

shown.   As  we  have  previously  discussed,  the  motion  implied 

must  have  a  short  coherence  length  and  presumably  gives  rise  to 

that  portion  of  x"  represented  by  (35).   The  remainder  of  the 

sum  in  (55')  is  not  small.   In  fact,  it  would  diverge  at  large 

( 2-dimensional )  k  if  not  cut  off  at  some  characteristic  length 
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Writing    k  =  (2 w/L) (K^ , K^)   fcr  integer   K^,  K^,       setting 
'^'rrax  "  (  2  Tr/a-Q )  ( ?  tt/L  )  =  L/Bq  ,    and  approximating  the  sum  in 
(55')  by  an  integral,  we  -have 


2    1 
X   =  X, 


(56) 


v;here  Xp   is  the  k  =  0  contribution.   Thus  the  "capillary  waves" 

k  /  0   make  3  large  (and  it  seems  ultimately  divergent)  contribution 

to  the  broadening  of  the  density  profile. 


5.   Surface  Dynamics 

We  have  been  led  to  a  picture  of  the  phase  transition 
region  in  which  a  cuasi-plane  symmetric  pattern  of  parallel 
surface  motions  of  relatively  short  coherence  length  is  further 
perturbed  by  waves  excited  in  the  effectively  elastic  surface. 
Our  analysis  has  been  partly  heuristic  and  partly  rigorous, 
raising  at  least  as  many  questions  as  it  claims  to  answer.   For 
example,  the  analysis  of  {^8)    applied  to  the  density  profile 
itself  has  been  seen  to  lead  to  a  divergence  going  as  the 
logarithm  of  the  system  size.   We  have  therefore  begun  a  detail- 
ed study  of  the  dynamics  of  the  surface  motions  in  computer 
simulation  . 

To  start  with,  we  have  plotted  in  Figures  11,  12,  and   13 
the  atomic  configurations  in  the  vapor,  surface,  and  liquid 
regions  at  one  particular  time.   Each  atom  is  represented  by  a 
circle  with  diameter  a  and  the  broken  circles  represent  the 
hidden  atoms  when  looking  from  the  vapor  region  into  the  surface, 
These  pictures  do  show  a  clustering  of  atoms  in  the  surface  with 
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well  defined  regions   of  hiph  and  low  density.   When  we  follow 
these  pictures  in  tirre,  we  observe  the  clusters  ipoving,  but  it  is 
difficult  with  our  data  to  distinguish  diffusive  from  wave-like 
propagation  of  the  clusters.   However,  over  as  few  as  800  con- 
figurations, the  irean  density  has  already  achieved  a  form  sensibly 
independent  of  x  and  y. 

To  examine  the  surface  motion  in  greater  detail,  it  is  first 
necessary  to  define  the  surface.   For  a  meaningful  microscopic 
visualization,  this  is  best  done  by  dividing  the  system  into  gas 
atoms  and  liouid  atoms.   For  a  rather  sharp  criterion,  we  may 
compare  the  local  enerpy  of  a  particle  with  some  reference  energy: 


r.  in  gas  if      E 
in  liouid  if 


I       «(r.-r.)  >  E 


(57) 


Indeed,  the  surface  obtained  from  the  liquid  surface  atom-atom 
connections  is  quite  insensitive  to  the  threshold  value  £„, 
reinforcing  the  picture  of  a  sharp  short  time  liquid-gas  inter- 
face.  This  is  shown  in  Figures  1 ^    and  15,  essentially  histograms 
of  C(x,y).   More  important  is  the  dynamical  information  in  the 
form  of  time-dependent  correlation  functions   F(k,t),   obtained 
by  replacing   r.  in  equation  C^O)    by  r.(t).   Here,  preliminary 
computations  at  least  show  that  F ( k  ,  t )   decays  to  equilibrium 
much  more  slowly  in  the  surface  repion  than  in  the  bulk,  evidence 
again  of  the  presence  of  slow  long  waves  manifested  presumably 
as  slow  large  clusters.  A    more  detailed  analysis  using  longer 
simulations  and  larger  surface  area  is  now  being  undertaken  and 
will  be  reported  in  due  course. 
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6.   Conclusions 
The  density  profile  fourd  by  Rao  and  Levesaue   appears 
superficially  to  indicate  a  sn^ooth  and  uninteresting  transition 
frorr  liquid  to  vapor  extending  over  a  broad  range  of  about 
five  length  units.   Put  we  find  as  a  result  of  the  further 
analysis  presented  here  that  the  interface  has  considerable 
structure.  In  particular,  the  transverse  structure  factor  does 
not  interpolate  sirioothly  in  the  interfacial  region.   On  the 
contrary,  it  exhibits  a  striking  peak  for  small  k  in  the  inter- 
face.  We  interpret  this  as  evidence  of  capillary  waves  in  the 
interface  whose  leading  terrr  for  srrall  k  gives  rise  to  a  struc- 
ture factor  contribution  S(z,z';k)  a  d'(2)p'(z').   In  addition, 
we  find  that  purely  longitudinal  correlations  are  consistent 
with  a  model  in  which  the  interface  is  taken  to  be  planar, 
sharp,  and  fluctuating  in  location.   Taken  with  the  evidence 
provided  by  snapshots  of  atoms  in  the  interface,  we  conclude 
that  the  delineation  between  liquid  and  vapor  is  rather  abrupt 
—  of  the  order  of  one  atomic  diameter  for  our  parameters  --  but 
that  this  boundary  fluctuates  markedly  in  position  and  time. 
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Captions  for  Figures 

1.  Density  Profile  for  Model  Argon  at  84°  K  (units  as  in  text). 

2.  Full  Density  Profile  for  84°  K  Model  Argon  in  Periodic 
Rectangular  Parallelepiped. 

3.  Typical  Sub-critical  Isotherm  for  Liquid-Vapor  Transition. 

4.  Radial  Distribution  Function   ggCr)  for  84°  K   Model 

Argon  at  Liquid  Density  n„  =  0.76. 
(2) 


Denote  az  =  1.38,   Open  Circles  Az  =  2.31,  Open  Triangles 
Az  =  3.23. 

6.  Density  Profile  from  Equation  (16*). 

7.  k-Dependence  of  Transverse  Structure  Factor  in  Liquid, 
z  =  0. 

8.  k-Dependence  of  Transverse  Structure  Factor  at  Mid-interface, 
z  =  7.5.   Curve  Shows  System  of  1728  Particles, 

Circles  Denote  System  of  1C24  Particles. 

9.  z-Dependence  of  Excess  Transverse  Structure  Factor  at 

k  =  2Tr/13-15,  AZ  =  0.4  (Line);     p(z)  is  shown  by  triangles. 
Comparison  Curve  of  I0p'(z)^  shown  by  Circles. 

10.  Low-k  Transverse  Correlations  at  Mid-interface. 
Solid  Curve  from  Equation  (47),  Computer  Simulation 
Points  Denoted  by  Crosses. 

11.  Typical  Configuration  of  Atoms  in  Vapor  Region,  Centers 
in  Slab  of  Thickness  az  =  1.0. 

12.  Typical  Configuration  of  Atoms  in  Interface  Region, 
Centers  in  Slab  of  Thickness   az  =  1.0. 

13.  Typical   Configuration  of  Atoms  in  Liquid  Region, 
Centers  in  Slab  of   Thickness  az  =  1.0. 

14.  Typical  Surface  Configuration  of  Atoms,  with  Energy 
Criterion  Eq  =  0. 

15.  Typical  Surface  Configuration  of  Atoms  with  Energy 
-.5. 
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Figure   6 
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Figure    7 
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Figure    8 
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Figure   9 
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Figure    10 
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Figure    11 
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Figure    12 
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Figure    13 
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